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My  research  works  under  the  support  of  the  ONR  grant  can  be  classified  into  two  major 
areas.  One  is  to  develop  an  adaptive  base  finite  element  method  for  computing  multiple  scale 
solutions  arising  from  material  science  and  turbulent  transport  [7,8].  The  other  is  to  develop 
stable  and  efficient  numerical  methods  to  compute  free  surface  flows  in  two  and  three  space 
dimensions  [1]. 

1,  An  adaptive  base  finite  element  method  for  multiple  scale  solutions 

Many  problems  of  fundamental  and  practical  importance  have  multiscale  solutions.  Com¬ 
posite  materials  and  turbulent  transport  in  high  Reynolds  number  flows  are  examples  of  this 
type.  Due  to  the  limited  computational  resources,  it  is  usually  impossible  to  completely 
resolve  all  small  scale  components  in  the  solution.  On  the  other  hand.  From  an  engineering 
perspective,  it  is  often  suflicient  to  predict  the  macroscopic  properties  of  these  systems,  such 
as  effective  conductivity,  elastic  moduli,  permeability,  and  eddy  diffusivity,  and  to  capture 
the  averaged  effect  of  small  scales  on  the  large  scales.  The  challenging  question  is:  Can  we 
capture  the  large  scale  property  correctly  without  completely  resolving  all  the  small  scale 
components? 

In  this  work  [7,8],  we  develop  a  new  method  of  computing  multi-scale  and  multicomponent 
problems  in  composite  materials  and  turbulent  transport  problems.  The  new  method  is  based 
on  a  general  theory  of  G-convergence  of  differential  operators  with  multiple  scale  structure. 
In  this  method,  the  base  functions  are  constructed  in  such  a  way  that  they  contain  the 
local  microscopic  properties  of  the  differential  operators.  Typically  the  size  of  the  element  is 
larger  than  a  certain  cut-off  scale  of  the  oscillatory  coefficient,  information  at  smaller  scales 
are  built  into  the  base  functions.  The  small  scale  information  in  the  base  functions  are 
brought  into  the  large  scale  (i.e.,  grid  scale  or  larger)  solution  through  the  coupling  of  the 
global  stiffness  matrix.  Thus,  the  large  scale  solution  is  correctly  captured.  In  this  sense,  we 
say  that  the  base  functions  are  adaptive  to  the  differential  operator. 

In  general,  the  adaptive  base  functions  need  to  be  constructed  numerically  except  for 
certain  special  cases,  such  as  piecewise  constant  coefficients  or  space-separable  problems. 
Since  the  construction  is  a  local  operation  within  the  elements,  it  can  be  done  in  perfectly 
parallel.  In  effect,  we  break  a  large  scale  computation  into  many  smaller  and  independent 
pieces;  this  is  in  contrast  to  traditional  domain  decomposition  methods,  in  which  the  pieces 
are  coupled  together.  Thus,  the  method  is  automatically  adapted  to  parallel  computers. 
Moreover,  once  the  small  scales  information  within  an  element  (a  base  function)  is  gathered 
into  the  global  stiffness  matrix,  the  computer  memory  used  for  those  base  functions  can 
be  reused  for  constructing  bases  of  the  next  element.  The  ability  to  gather  the  fine  scale 
information  in  such  a  sequential  manner  leads  to  drastic  saving  of  computer  memory.  This 
makes  it  possible  for  the  method  to  tackle  much  larger  practical  problems  than  the  direct 
numerical  methods.  Clearly,  when  the  small  scale  features  are  localized  in  a  region,  we  need 
only  to  construct  the  adaptive  bases  in  that  region.  Therefore,  the  adaptive-base  method  can 
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be  used  as  an  alternative  to  the  h-p  finite  element  method  or  the  adaptive  mesh  refinement 
methods  when  the  averaged  effect  of  the  small  scales  are  of  the  main  concern.  Another  feature 
of  our  approach  is  that  it  can  handle  general  multi-scale  problems  without  the  requirement  of 
scale  separation,  a  property  which  is  very  important  for  practical  applications.  Our  numerical 
experiments  convincingly  demonstrate  this  effect. 

Convergence  analysis  of  our  adaptive  finite  element  method  is  a  difficult  task.  Analysis 
for  conventional  finite  element  methods  can  not  be  used  to  analyze  our  adaptive  method.  A 
more  refined  analysis  which  takes  into  account  the  microscale  structures  of  the  solution  is 
required  to  reveal  the  convergence  property  of  our  method.  In  the  case  of  periodic  multi-scale 
oscillatory  coefficients,  we  have  been  able  to  demonstrate  the  convergence  of  our  method.  Our 
analysis  shows  that  the  adaptive  method  converges  independent  of  the  small  scale  feature  of 
the  solution  [8] .  In  a  subsequent  paper,  we  have  carried  out  a  number  of  careful  performance 
studies  using  massively  parallel  super  computers.  Comparison  with  other  existing  numerical 
methods  is  also  made.  Our  study  shows  that  our  method  is  very  competitive  with  other 
numerical  methods.  Our  method  can  compute  multiscale  solutions  of  small  scales  that  are  one 
or  two  orders  magnitude  smaller  than  scales  that  can  be  computed  by  conventional  methods. 
And  the  overhead  in  computing  the  oscillatory  bases  function  only  double  the  amount  of 
works  required  for  a  conventional  finite  element  method.  To  demonstrate  the  adaptivity 
and  the  robustness  of  our  method,  we  have  computed  random  oscillatory  coefficient  with 
fractal  dimensions  for  practical  problems  in  composite  materials  and  porous  media  flows. 
Convergence  is  observed  in  all  our  numerical  experiments  [7].  In  the  future,  we  plan  to  use 
our  method  to  study  the  effective  conductivity  of  certain  composite  materials  with  large 
aspect  ratio  and  many  unseparable  small  scale.  We  will  also  study  wave  propagation  in 
random  media  and  the  effective  diffusivity  and  possible  anomalous  diffusion  in  turbulent 
transport  problems. 

2.  Efficient  and  stable  numerical  methods  for  free  surface  flows 

Many  physically  interesting  problems  involve  propagation  of  free  surfaces.  Thin  film 
growth  for  manufacturing  materials  for  electronic  devices,  solidification  problem  for  crystal 
growth,  Hele-Shaw  cells  for  pattern  formation  are  some  of  the  significant  examples.  These 
problems  present  a  great  challenge  to  physicists  and  applied  mathematicians  because  the 
underlying  problem  is  very  singular.  The  physical  solution  is  sensitive  to  small  perturbations. 
Straight-forward  discretizations  may  lead  to  numerical  instabilities.  Moreover,  in  many 
situations,  the  numerical  modeling  is  especially  difficult  due  to  the  presence  of  topological 
changes  or  formation  of  singularities  in  the  surface.  For  example,  in  crystal  growth  and 
thin  film  growth,  an  initial  smooth  front  can  develop  cusps  and  crack-like  singularities,  and 
isolated  islands  of  film  material  can  merge.  Another  difficulty  is  due  to  the  presence  of 
surface  tension  in  the  free  surface.  It  introduces  a  very  stiff  term  through  the  local  mean 
curvature.  In  the  case  of  thin  film  growth,  the  situation  is  even  worse.  The  normal  velocity 
of  the  free  surface  is  proportional  to  the  second  derivative  of  the  local  curvature.  This  would 
pose  an  extremely  severe  stability  constraint  for  time  integration  if  an  explicit  scheme  is 
used. 

In  collaboration  with  Lowengrub  and  Shelley  [2],  we  introduce  a  novel  idea  to  remove 
the  stiffness  of  surface  tension  by  reformulating  the  problem  in  the  tangent  angle  and  the 
arclength  metric  variables.  Using  this  reformulation,  we  can  decompose  the  small  scale 
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features  and  integrate  the  most  singular  part  of  the  equations  implicitly.  Our  reformulation 
enables  us  to  invert  the  implicit  solution  explicitly  by  FFT,  this  greatly  improves  the  speed 
of  the  time  integration.  For  computations  of  the  vortex  sheet  roll-up  in  an  Euler  flow  with 
surface  tension  using  a  modest  number  of  points  (128),  the  time-step  can  be  chosen  250  times 
larger  than  for  an  analogous  explicit  method.  Many  problems  which  were  previously  un¬ 
obtainable  now  become  tractable  using  our  techniques  and  new  phenomena  were  discovered. 
Using  our  reformulated  method,  a  new  type  of  topological  singularity  was  found  numerically 
for  fluid  interfaces  with  surface  tension,  which  is  different  from  the  curvature  singularity  due 
to  the  Kelvin-Helmholtz  instability.  A  careful  numerical  study  has  been  recently  carried 
by  us  to  study  the  scaling  property  of  this  pinching  singularity  using  a  high  order  adaptive 
version  of  our  reformulated  boundary  integral  method  [3].  A  number  of  interesting  physical 
phenomena  are  revealed.  Numerical  study  of  the  two-density  unstably  stratified  flow  with 
surface  tension  has  also  been  carried  in  a  joint  work  with  my  student,  Hector  Ceniceros  [4]. 
Our  study  shows  that  a  similar  pinching  singularity  is  also  observed  in  the  roll-up  of  the  free 
surface.  The  type  of  singularity  is  qualitatively  similar  to  the  ones  we  observed  for  the  vortex 
sheet  with  surface  tension.  Convergence  of  the  reformulated  boundary  integral  method  has 
been  obtained  using  delicate  energy  estimates. 

Another  significant  result  we  obtained  is  to  develop  and  analyze  a  stable,  high  order  3-D 
boundary  integral  method  for  computing  water  waves  [9,10].  The  stability  analysis  of  a  3-D 
boundary  integral  method  is  much  more  difficult  and  subtle  than  its  2-D  counterpart.  For 
a  3-D  problem,  the  singular  kernel  has  an  un-removable  branch  point  singularity.  Moreover, 
it  is  much  more  difficult  to  obtain  the  spectral  property  of  the  singular  integral  operators 
with  variable  coefficients.  In  our  analysis,  we  develop  a  framework  under  which  spectral 
properties  of  various  singular  operators  can  be  analyzed  precisely.  We  found  that  the  3-D 
boundary  integral  method  is  generically  unstable.  Using  a  delicate  stability  analysis,  we 
show  how  to  remove  this  numerical  instability  and  construct  a  stable  method  for  computing 
3-D  water  waves. 

In  an  effort  to  compute  through  topological  singularity  beyond  the  pinching  of  the  inter¬ 
face,  we  derive  a  level  set  formulation  for  incompressible,  immiscible  Navier-Stokes  equations 
separated  by  a  free  surface  [5].  This  is  done  in  collaboration  with  Chang,  Merriman  and 
Osher.  The  flow  we  consider  has  discontinuous  density  and  viscosity.  The  effect  of  surface 
tension  is  also  included.  Based  on  this  formulation,  a  second  order  projection  method  can  be 
used  to  approximate  the  evolution  equations.  This  approach  can  be  considered  as  a  method 
of  front  capturing  type  since  no  explicit  information  about  the  free  surfaces  is  required  in 
the  solution  procedure.  The  free  surface  is  recovered  at  the  end  of  the  computation  by  lo¬ 
cating  the  zero  level  set  of  a  smooth  function.  This  numerical  method  is  efficient  and  is 
capable  of  simulating  incompressible  flow  where  change  of  topology  in  the  fluid  interface 
occurs,  such  as  merging  and  reconnection.  An  area-preserving  re-initialization  procedure  is 
introduced  to  prevent  loss  of  mass  due  to  numerical  diffusion.  This  property  is  essential  for 
many  practical  applications.  Recently,  we  (with  Osher)  have  obtained  an  improved  level  set 
formulation  for  vortex  sheets  with  surface  tension  [6].  In  this  formulation,  the  singular  delta 
function  is  factored  out  from  the  evolution  equation.  This  gives  a  more  robust  formulation 
for  both  theoretical  and  computational  study.  We  intend  to  perform  a  careful  numerical 
study  and  compare  with  the  result  obtained  by  using  our  reformulated  boundary  integral 
method.  The  advantage  of  using  the  level  set  method  is  that  we  can  compute  beyond  the 
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pinching  singularity  without  any  difficulty. 
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